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A Monte Carlo method based on a density-of-states sampling is proposed for study of arbitrary 
statistical mechanical ensembles in a continuum. A random walk in the two-dimensional space of 
particle number and energy is used to estimate the density of states of the system; this density 
of states is continuously updated as the random walk visits individual states. The validity and 
usefulness of the method are demonstrated by applying it to the simulation of a Lennard- Jones 
fluid. Results for its thermodynamic properties, including the vapor-liquid phase coexistence curve, 
are shown to be in good agreement with high-accuracy literature data. 
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INTRODUCTION 

The free energy landscapes of complex systems, such as 
proteins or glasses, are characterized by the existence of 
deep, local free energy minima. These minima pose sig- 
nificant obstacles for molecular simulations; once a sys- 
tem is trapped in a minimum, conventional algorithms 
are unable to explore configurations pertaining to other, 
relevant regions of phase space. Several Monte Carlo 
methodologies have been developed in the last decade 
to circumvent the sampling problem. Examples include 
expanded ensembles |], ||], multicanonical algorithms ||, 
and parallel tempering formalisms Gj. All of these tech- 
niques have improved considerably our ability to simu- 
late the equilibrium properties and structure of complex 
materials, including polymers, proteins in solution, or or- 
ganic glasses. 

Most of these techniques have relied on Metropolis' et 
al. original prescription^, in which trial configurations 
of a system are accepted or rejected according to proba- 
bility distributions pertaining to conventional statistical 
mechanical ensembles (or minor alterations thereof). A 
notable exception is provided by multicanonical meth- 
ods, in which trial configurations are accepted according 
to multicanonical weights constructed in such a way as 
to "flaten out" high energy barriers between distinct con- 
figurations. Unfortunately, the required multicanonical 
weights must be determined through a tedious and highly 
computationally demanding iterative process. In multi- 
canonical methods as well as in more conventional Monte 
Carlo techniques for molecular simulation, the condition 
of detailed balance is satisfied by construction. 

Recently, however, Wang and LandauQ proposed a 
promising approach for lattice systems which eliminates 
some of the shortcomings of the original multicanonical 
prescription. The key to that approach is to accept trial 
configurations of the system according to a running es- 
timate of the density of states; by design, their formal- 
ism does not strictly enforce detailed balance. A random 



walk in energy space is used to evenly visit each energy 
level. The density of states of an energy level is modi- 
fied by an arbitrary convergence factor when that energy 
level is visited. By controlling that factor in a systematic 
manner, these authors were able to generate the density 
of states of an Ising lattice system to high accuracy in a 
self-consistent way. We refer to this method as Density- 
Of-States (DOS) Monte Carlo. 

For simulations on a lattice, DOS Monte Carlo 
promises to offer significant advantages over previously 
available techniques. It is therefore of interest to explore 
whether analogous ideas can be used in a continuum, in 
the context of a realistic fluid. Several challenges must 
be overcome: first, random moves on a lattice system can 
only give rise to a small set of discrete energy changes, 
thereby simplifying considerably the nature of the simu- 
lations. Random displacements in a continuum result in 
unpredictable energy changes, and it is unclear whether 
DOS Monte Carlo can be implemented at all. Second, 
and perhaps more important, DOS Monte Carlo gener- 
ates an estimate of the density of states to within a mul- 
tiplicative constant. Such a constant depends on the par- 
ticular length of a simulation, on system size (volume), 
and on density (number of particles). On a lattice, it is 
common practice to study systems of constant composi- 
tion. Furthermore, the system is generally assumed to be 
incompressible. Real fluids, however, are not incompress- 
ible, and their study over any reasonable range of density 
would require that the absolute density of states (and its 
multiplicative constant) be known. A naive extension of 
the lattice DOS method to a real fluid would therefore re- 
quire that multiple simulations be conducted at different 
densities, and that the resulting densities of states corre- 
sponding to each calculation be matched with each other 
in some fashion to estimate each multiplicative constant. 
This procedure would be prone to uncertainties related 
to the manner in which different thermodynamic states 
are combined, it would be time consuming, and it would 
not offer significant advantages over thermodynamic in- 



2 



tegration or histogram reweightingQ. 

In this work we propose a DOS formulation which ad- 
dresses these issues. Its implementation is demonstrated 
in the context of a Lennard- Jones fluid, for which high- 
accuracy thermodynamic data area available. It is shown 
that the proposed simulation technique is able to gener- 
ate the partition function of the system over a wide range 
of energy and density, which comprises an infinitely di- 
lute gas and a dense liquid. This is achieved within a sin- 
gle simulation, at a small fraction of the computational 
demands of currently available simulation techniques. 



DENSITY OF STATES MONTE CARLO 

On a lattice, a DOS Monte Carlo method generates a 
random walk in energy space by flipping spins in a ran- 
dom manner. A trial flip is accepted according to crite- 
ria which eventually result in a flat distribution of energy. 
Generating a flat energy distribution requires that micro- 
scopic states having internal energy E be visited with a 
probability inversely proportional to the density of states, 
g{E). Accepting trial spin flips with probability 



p(Ei — > E2) — min 



g(E 2 ) 



(1) 



(where E\ and E2 are the energies before and after the 
spin is flipped), would therefore result in a flat energy 
distribution. 

Unfortunately, the density of states g(E) is not known 
a priori. The central idea in DOS Monte Carlo is to 
construct the density of states on the fly. At the begin- 
ning of a simulation, a constant density of states is as- 
sumed for all energy levels. As the simulation proceeds, 
the instantaneous values of g(E) are used to accept trial 
flips. Every time that an energy level E is visited, the 
density of states for that level is modified according to 
g(E) — > g(E)* /, where / > 1 is an arbitrary convergence 
factor. An energy histogram H(E) is accumulated during 
the course of the simulation. When that histogram be- 
comes sufficiently flat, the convergence factor is reduced 
to a finer value, e.g. = \ffi\ the histogram H(E) is 
reset to zero, and the simulation is continued. This pro- 
cess is repeated until the convergence factor / becomes 
arbitrarily small (say, less than exp(10~ 8 )). 

In the many-body fluid formulation proposed here, the 
number of particles or molecules in the system is allowed 
to fluctuate. We consider a phase space characterized by 
internal energy E and number of particles N . To sample 
relevant regions of phase space efficiently, our strategy is 
to have each pair of (N, E) points be visited uniformly. 
In other words, a microscopic state having N particles 
and energy level E should be visited with probability 
l/g(N, E), where g(N, E) represents the two-dimensional 
degeneracy of (N, E) states. 



For simplicity, only two types of trial moves are con- 
sidered. The first type of move consists of simple trial 
displacements, in which the coordinates of a randomly 
chosen particle are altered by a small random amount. 
A trial move is then accepted or rejected according to 



p{E 1 ,N -> E 2 ,N) = min 



gjEuN) 
g(E 2 ,N) 



(2) 



The second type of move consists of trial insertions or 
destruction of particles. For a trial insertion, a particle 
is introduced into the system at a random position. This 
move is accepted with probability 



p(N -> N+l) = min 



1, 



V 



g(N,E° ld ) 



{N + l)A 3 g(N+l,E nev/ ) 



(3) 

where V is the volume of the system, and A is the de 
Broglie thermal wave length. In Equation (§), E old and 
E ncw represent the energies of the system before and af- 
ter a particle is introduced, respectively, and g(N,E old ) 
and g(N + l,E ncw ) are the corresponding densities of 
states. For a trial destruction, a randomly chosen parti- 
cle is removed from the system, and the move is accepted 
with probability: 



p(N ^ N -1) 



1. 



g(N,E° 



Id \ 



V 



g{N - l,i; new ) 



(4) 



It is of interest to remark that traditional, Metropolis- 
type simulation techniques penalize trial insertions or de- 
structions of particles by a factor proportional to the ex- 
ponential of the energy, exp(— (3AE), where AE repre- 
sents the change in energy created by the trial move, 
and where (3 = X/k^T (T is the temperature). If AE 
exceeds a few k^T, the trial move is generally rejected. 
That factor is absent from the algorithm proposed in this 
work, thereby facilitating considerably the creation or de- 
struction of molecules and the sampling of high-density 
configurations. 

A running estimate of the two-dimensional density of 
states is continuously updated as the simulation pro- 
ceeds. When a configuration having N particles and 
energy E is visited, the current value of g(N,E) is mul- 
tiplied by a convergence factor /. A two-dimensional 
histogram of number of particles and energy H(N,E) is 
constructed. When that histogram is deemed to be suffi- 
ciently flat, the histogram is discarded and the simulation 
is continued, this time with a smaller convergence factor. 

Having generated a density of states according to the 
procedure outlined above, a partition function for the 
system can be constructed at any given temperature and 
chemical potential according to 



HCr» = ]T]r. g (A,£;) e 



(5) 



N E 



The partition function determined from Equation ^| is 
only known to within an arbitrary constant multiplier. 
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Our simulation, however, comprises the special case N = 
0; the density of states for that case is known to be unity, 
thereby providing a means to the determine the absolute 
value of 3 for all other states. Thermodynamic prop- 
erties of interest can subsequently be determined from 
knowledge of 5(T, fx). For example, the thermodynamic 
pressure and the thermodynamic internal energy can be 
calculated according to 



P (T>) = ^io g s(T, M ) 



(6) 



and 



£/(!» = 



N E 



-llE+NPv 



(7) 



JV E 



The compressibility of the fluid can also be inferred 
from g(N,E) by considering fluctuations in the number 
of particles according to 



V (N 2 ) - (N)' 



K{T,fi) v{dpJ Ai k B T iXr 



(8) 



where 



■PE+NPti 



(N) 



N E 



and 



(A 2 ) = 



N E 



N E 



(9) 



-pE+NPn 



(10) 



N E 



In many applications, it is of particular interest to de- 
termine the precise location of first or second order ther- 
modynamic phase transitions. In the particular case of a 
simple fluid, to calculate the liquid-vapor binodal curve, 
a simple two-state construction can be used. At a tem- 
perature well below the critical point, the equilibrium 
density distribution exhibits two distinct peaks; a thresh- 
old number of particles Nq can be designated, such that 
all N < No states can be regarded as pertaining to a 
"vapor" branch, and all N > No states as belonging to 
a "liquid" branch. The pressures corresponding to these 
two branches can be calculated according to 



k B T 
V 

k B T 
V 



l0 § E ^g{N,E)e-^+^" (11) 

N<N E 

l0 S E T,9(N,E)e^ E+N ^. (12) 

JV>Af E 



For any given temperature, a phase coexistence point 
can be found by carefully tuning the value of chemical 
potential /io in such a way as to satisfy the condition 
p v (T,p ) = P L {T,fio)- The coexistence densities of the 
vapor and liquid phases can then be calculated as 



(3E+N[3^o 



p v (T) 



N<N E 



v E E^ £ ) e " 



■PE+NPho 



N<N E 



J2 Y,Ng(N,E)e 



-pE+NPfio 



p L {T) — N>No E 



N>N E 



(13) 



(14) 



Note that the equilibrium pressure obtained from 
Equations ll| and |l2| differs from that obtained from 
Equation p by a term — k B T log2/V, which arises from 
finite-size effects. For systems having a large enough vol- 
ume, however, this difference is negligible. 

One particular problem that must be addressed in the 
method proposed above is that the relevant range of en- 
ergy is strongly dependent on the number of particles 
in the system. The energy levels accessible to a two- 
particle system are different than those accessible to a 
100-particle system. To determine energy ranges for dif- 
ferent system sizes, two short preliminary DOS simula- 
tions are run, one at the lowest temperature, and the 
other at the highest temperature. The goal of these simu- 
lations is to achieve a flat distribution of N, as opposed to 
trying to make the (N, E) histogram flat; the energy dis- 
tribution is dictated by a conventional Boltzmann weight. 
A flat N distribution would be obtained if each micro- 
scopic state with N particles was visited with probability 
l/Q(N), where Q(N) is the canonical partition function 
(of a system having N particles) . 

The scheme followed in these preliminary runs is simi- 
lar to that adopted above: a table of partition functions 
is constructed, with entries for each value of N. At the 
beginning of a simulation, the partition functions are set 
to unity for all entries of N. Two types of moves are 
used: particle displacements, accepted according to con- 
ventional Metropolis criteria, and particle insertion or 
destruction moves, accepted with probability 



p(N -> N+l) = min 
and 

p(N -> N-l) = min 



1, 



V 



« „ Q ( N \, x exp(-/3A£) 



(Ar + 1)A 3 Q(N + 1) 
Q(N) 



(15) 

x exp(-(3AE) 
(16) 

where AE is the energy change associated with the trial 
particle insertion or destruction. Upon each visit to a 
microscopic state, the corresponding Q(N) is updated by 



iVA 3 

1, x 

' V 



Q(N-l) 



4 



multiplying it by an arbitrary convergence factor /. The 
minimum and maximum values of energy corresponding 
to each number of particles are tracked during the sim- 
ulation; the minimum energy at the lowest temperature 
and the maximum energy of the highest temperature de- 
fine the relevant energy range to be used in the two- 
dimensional g(N,E) DOS production simulations. 

Since the purpose of the preliminary runs is to de- 
termine the energy range, it is not necessary to gener- 
ate a perfectly flat histogram. The only requirement is 
that each number of particles be visited with enough fre- 
quency. It is important to emphasize, however, that the 
scheme described above for Q(N) is very useful on its own 
right, for example for calculation of the chemical poten- 
tial of a fluid. If this DOS Q(N) simulation is run until 
convergence is achieved (i.e. until the N histograms are 
flat and the convergence factor / is small enough), the 
result is the free energy of the system as a function of 
N. This information is particularly useful in expanded 
ensemble simulations 0, 0, where the weights associated 
with individual expanded states (which are closely re- 
lated to free energies) must be determined before a pro- 
duction simulation. In the scheme described here, there 
is no need to determine those weights, as they would be 
determined directly through the course of the simulation. 

APPLICATION TO LENNARD-JONES FLUIDS 



an energy bin size of e is used to construct the density 
of states and the required histograms. Approximately 
5 x 10 7 Monte Carlo steps were used to generate the 
complete density of states. 

Figure Q shows the density of states as a function of en- 
ergy and number of particles. Note that for most (N, E) 
pairs, the corresponding density of states is less than 
unity. This apparently paradoxical feature is actually 
due to the fact that the de Broglie thermal wave length 
is set to unity in our simulations. 

Figure || shows the phase diagram of the truncated 
Lennard-Jones fluid. The solid line was calculated from 
the density of states determined in this work; the trian- 
gles show literature results for the same system || . The 
agreement between the two sets of data is good. The fig- 
ure shows several isobars, also calculated from the den- 
sity of states. These isobars are in excellent agreement 
with pressure calculations from conventional canonical- 
ensemble simulations (results not shown). The isother- 
mal compressibility of the system, calculated according 
to Equation ||, is shown in Figure || as a function of den- 
sity along several isobars. These results are also consis- 
tent with those of conventional simulations. A distinct 
peak can be observed in the compressibility near p = 0.3. 
As expected, the peak becomes more pronounced as the 
pressure approaches the critical pressure of the system 
(and fluctuations become more prominent). 



The remainder of this letter describes results from the 
application of the DOS method outlined above to the 
particular case of a truncated Lennard-Jones fluid. This 
model exhibits most of the main features that one expects 
to find in most realistic fluids, and offers the advantage 
that high-accuracy simulation data are available for its 
thermodynamic properties. The potential energy of in- 
teraction between two particles is of the form 

C r>r c 

U{r) ^\M& 2 -&} ^> 

where r is the distance between the particles, and r c is the 
cutoff distance. To compare our results to those reported 
in the literature, we use r c = 2.5a. The box length of 
the system is set to L = 5er. 

Before a simulation is conducted, a relevant range of 
interest for the number of particles and the energy must 
be specified. In this work, that range is set to be between 
to 110 for the number of particles, which covers the den- 
sity range p* — 0-0.88. The range of energy comprises 
E/e = —690 at one end of the spectrum, and E/e = 10 
at the other end; this range corresponds to temperatures 
between 0.5 < T* < 1.5. 

In contrast to a spin lattice, the potential energy of 
a Lennard-Jones fluid is continuous; a discretization of 
the energy must therefore be introduced. In this work, 



CONCLUSIONS 

A density of states Monte Carlo method has been pre- 
sented for simulation of the thermodynamic properties 
of realistic fluids over wide ranges of density and en- 
ergy. This method permits calculation of virtually all of 
the thermodynamic properties of a system, including its 
phase behavior, as a function of density or temperature; 
all of this information is generated from a single simu- 
lation. Furthermore, by virtue of the way in which trial 
configurations are generated, states of the system that 
would otherwise be difficult to sample (e.g. low temper- 
atures or high densities) can be studied efficiently and 
with remarkable accuracy. An additional benefit of the 
method is that prior knowledge of the behavior of the 
fluid (such as approximate location of phase boundaries) 
is not required; a "blind" simulation is able to sample 
relevant configurations of the system with little, if any 
user-supplied guidance. 

For the particular case of a truncated Lennard-Jones 
fluid, it has been shown that the procedure proposed in 
this work generates thermodynamic information in quan- 
titative agreement with high accuracy literature data for 
the same fluid. The coexistence curve, vapor pressure, 
and isothermal compressibility are all in agreement with 
results of conventional simulation techniques. For this 
fluid, the entire range of temperature and density con- 
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sidered in this work can be generated in several hours of 
computer time. 

While in this work the DOS Monte Carlo method has 
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FIG. 1: Two-dimensional density of states of the truncated 
Lennard- Jones fluid. Different lines correspond to a differ- 
ent number of particles. The number of particles increases 
monotonically from right to left. 
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FIG. 2: Phase diagram of the truncated Lennard- Jones fluid. 
The solid line shows the results of this work; the triangles 
depict literature data for the same system The dashed 
lines are isobars calculated from the density of states. 




FIG. 3: Compressibility of the truncated Lennard- Jones fluid 
as a function of density, along various supercritical isobars. 



